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Abstract 

The investigation of freezing transitions of single polymers is computationally demanding, since surface effects dominate 
,the nucleation process. In recent studies we have systematically shown that the freezing properties of flexible, elastic 
polymers depend on the precise chain length. Performing multicanonical Monte Carlo simulations, we faced several 
computational challenges in connection with liquid-solid and solid-solid transitions. For this reason, we developed novel 
methods and update strategies to overcome the arising problems. We introduce novel Monte Carlo moves and two 
extensions to the multicanonical method. 
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1. Introduction 

Induced by the rapidly increasing efficiency and avail- 
ability of computational resources, the field of compu- 
tational physics has gained tremendously in importance 
within the last decades, and it is today regarded as physics' 
third pillar alongside experimental and theoretical physics. 
In addition to the innovations in hardware, simulation 
techniques have evolved further, and in fact, the greater 
improvements have resulted from better methods rather 
than from faster computers. A particularly important ap- 
plication is the investigation of thermodynamic properties 
of complex systems by means of Markov chain Monte Carlo 
methods. Starting sixty years ago with the Metropolis al- 
gorithm [l[, which emulates the canonical ensemble, the 
arsenal of algorithms has been extended and more sophis- 
ticated methods have been introduced. Among the most 
powerful simulation techniques are generalized-ensemble 
methods such as parallel tempering 0, Q , multicanonical 
samplinglj , simulated tempering [5[ , or the Wang-Landau 
method [6J], which allow in principle to collect all infor- 
mation about the entire thermodynamic behavior of the 
investigated system in a single simulation. However, de- 
pending on the considered system, substantial difficulties 
can occur, part of which are specifically related to prop- 
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erties of the system being studied, whereas others, like 
broken ergodicity, are of more general nature. 

In a recent study on flexible homopolymers 0, H| , we 
encountered a number of problems of both kinds and devel- 
oped new simulation techniques to overcome these. Some 
of them are rather specific to polymers, while others are 
more general and can also be applied to nonmolecular sys- 
tems. Combining our strategics we were able to boost the 
efficiency of our algorithms and to perform very precise 
simulations of systems which could not be investigated in 
this quality before. 

The purpose of this paper is to explain our methods in 
detail. After a short introduction of the applied polymer 
model in the next section, we briefly explain in section 
3 the multicanonical Monte Carlo method, which served 
as the basic algorithm in our simulations. The following 
section 4 is dedicated to the applied conformational up- 
dates and includes a new general optimization strategy for 
basic updates of systems with continuous degrees of free- 
dom. Afterwards we introduce and motivate in section 5 
two general extensions to the multicanonical method, and 
finish in section 6 with some concluding remarks. 

2. Model 

In our simulations we employed a bead-spring model 
for flexible, elastic polymers. For a specified set of monomer 
coordinates {X}, the energy of a polymer conformation is 
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given by 



for which the most common solution is given by 



JV-l JV 

E({X})=J2 E ^(IXi-Xyl) 
i— 1 

N-l 

+ ]TP b (|X m -X l |). 

i=l 

Here, the non-bonded interaction 

-Enb(r) = £ , Lj(min{r,r c }) - £ , L j(r c ) 



(1) 



(2) 



corresponds to a truncated and shifted Lcnnard- Jones (LJ) 
potential 

P LJ (r)=4[(a/r) 12 -(a/r) 6 ] (3) 
Pairs of bonded monomers fur- 



with the cutoff radius r, 
ther interact via 

K 



Mr) 



R 2 ln(l-[(r-r )/R] 2 ), 



(4) 



which is the standard finitely extensible non-linear clastic 
(FENE) potential. The parameters are chosen such that 
the minima of both potentials coincide at ro, in order to 
prevent frustration. For details of the paramctrization see 

fig. 

This model belongs to the class of coarse-grained mod- 
els, i.e., microscopic details have been traded for generality 
and handiness. However, accurate simulations are still a 
substantial challenge. 

3. Multicanonical Monte Carlo Sampling 

Before we discuss our novel simulation strategies, let us 
first recall basic principles of Markov chain Monte Carlo 
simulations (lpj , for which acceptance criteria are obtained 
from the master equation: 



/'->''j • 



(5) 



where P M (t) denotes the probability for a state li to occur 
at time t and W v ^^ is the transition probability from state 
v to li. In stationary equilibrium, where dP^(t)/dt = 0, 
this equation is solved by: 



II rv fJ.—H' 1 



(6) 



called "detailed balance" . The transition probability W y ^^ 
is the product of the probability of selecting the update 
proposal and the probability of accepting it: 

W v ^ = W^W^. (7) 

Symmetric selection probabilities 
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However, for convenience or increased sampling efficiency, 
it is useful to introduce Monte Carlo updates where the 
selection probabilities are unequal: 
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in which case 
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Then, the more general expression 
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of the acceptance probability is required. It has been 
demonstrated that such weighted updates can enable a 
much more efficient sampling of the system conformations 
llfl . compared with symmetrically chosen selection prob- 
abilities. This also applies to simulations in the grand- 
canonical ensemble (constant chemical potential) or a con- 
stant pressure in the Npt ensemble, where volume fluctu- 
ations are relevant 

The goal of the multicanonical method Q is to gener- 
ate a flat histogram H over a certain macroscopic observ- 
able which in our case is the energy E. This is achieved 
by introducing a weight function oj(E) which is inversely 
proportional to the density of states g(E): 



H(E) = const = w(E)g(E), 
u(E) cx g'\E). 



(14) 
(15) 



A single point in state space (conformation) fx = {X} is in 
the multicanonical ensemble represented by a probability 
density which is proportional to the weight function and 
is therefore depending only on the energy: 



P {x} cx W (P({X})). 



(16) 



The acceptance probability for a proposed Monte Carlo 
move is according to (flU|) 



IT' 



{X}^{X'} 



min 1, 



u(E({X'}))W. 



{X'}^{X} 



(17) 
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Usually, the density of states and hence the weight function 
is not known in the beginning and has to be estimated by 
iterative procedures such as error weighted accumulation 
fl3j ] or the Wang-Landau method [1] . 

4. Conformational Update Proposals 

4--1- Displacement move with energy dependent maximal 
step length 

When investigating many-particle systems by means of 
Monte Carlo simulations, the simplest possible conforma- 
tional update is the displacement of a single particle to a 



uniformly distributed random position within a sphere^ 
around its original location X$: 



X ■ — X, 



r, with Irj < r n 



(18) 



In the case of a flexible polymer with elastic bonds, such 
updates can be applied to all monomers. Thereby, the size 
''max of the sphere crucially influences the performance of 
the simulation. A larger sphere allows the system to per- 
form extended steps in conformational space and is there- 
fore appropriate for simulations at high temperatures. If 
the temperature is lowered, the efficiency decreases since 
the proposed steps arc now too large, and the system 
will not smoothly descend to narrow local energy min- 
ima. Moreover, if the system eventually finds an energy 
minimum, further moves are unlikely to be accepted, since 
the proposed changes will almost certainly result in a huge 
increase in energy. In consequence, smaller spheres should 
be used when a system with a rough energy landscape 
is investigated at low temperatures. It is simple to in- 
corporate variable sphere radii into simulation techniques 
such as Metropolis [l[ , parallel tempering Q , or simulated 
tempering [H[ by assigning suitable sphere radii to each 
temperature, i.e., to use 7" max (T) instead of r max , since 
for each of these methods a (sub) ensemble is associated to 
each single temperature and detailed balance is satisfied. 
Changes in temperature are usually performed separately 
from moves in conformational space and hence need not 
to be considered here. 

The situation is more complicated for multicanonical 
and Wang-Landau sampling, where a simulation temper- 
ature does not exist. Instead, the entire state space is 
sampled in a single generalized ensemble, making it dif- 
ficult to choose a single sphere radius that leads to ade- 
quate performance. However, the application of variable 
sphere radii is highly desirable, as it would greatly im- 
prove simulation efficiency. Since we require large steps at 
high and small steps at low energies, the energy itself ap- 
pears to be a well-suited control parameter for the sphere 
radii. However, using the standard multicanonical method 
with a maximal step length that depends on energy, and 
therefore changes in time, would cause a violation of the 
detailed balance condition. 

Let us discuss this in more detail by considering a 
displacement of the fcth monomer. Assume a conforma- 
tion {X' 1 } with a certain relatively high energy E^, and 
assume further, the maximum step length r max (Eh), is 
comparatively large. During the following update the sys- 
tem might jump to a rather small energy Ei with a much 
smaller sphere radius r max (Ei) < r max {Eh)- That means 
the maximum step length for the next update is smaller 
than for the first. As one consequence, the system some- 
times cannot reach the starting point Xjj within a single 
step, hence detailed balance is clearly violated. This is 



instead of a sphere, any three-dimensional body which is invari- 
ant under inversion of coordinates, e.g., an adequately oriented cube, 
would serve as well. 



the case if the distance between the two positions exceeds 



the smaller sphere radius |X£ — Xj 



.{Ei). Note 



XjJ < r mayi {Eh) holds by definition. Even if 



that |X£ 

this is not the case and the starting point lies within the 
smaller sphere, detailed balance is not fulfilled, because 
the probability densities for selecting the forward and the 
backward update are different and © is violated. Fortu- 
nately, according to (fT5|) . the emerging bias can easily be 
corrected. The probability density of proposing a certain 
displacement equals the inverse volume of the sphere: 



V( 

0, else 
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(19) 



For |X£ — Xl| < r max (E v ), one obtains according to (fT2")l 
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Hence, the final acceptance criterion reads 



min ( 1 



0, else. 



max I J - / l* 



(20) 



(21) 



Remember that the case |X£ — X£| > r max (E v ) cannot 
occur and is therefore not considered. 

In principle, any strictly positive function r max (E) can 
be employed, but here we are searching for a function 
that results in appropriate acceptance rates for all ener- 
gies. Therefore we start with a flat function and perform a 
tuning procedure. First, we apply a standard binning, i.e., 
we divide the energy axis in intervals in which r max (£') is 
constant, i.e., \{Ei<E< E i+ \ then r max (E) = r max (Ei), 
with a fixed interval size AE = -Ej+i — E^. The value of 
fma.x(Ei) shall now be adjusted such that approximately 
two third of all proposed updates increase the energy while 
the remaining third leads to lower energies. It is reason- 
able to assume that such values for r max (Ei) exist, since 
for very small values the accessible part of the energy land- 
scape resembles a tilted hypcrplane with one half belong- 
ing to higher and the other half to lower energies. If on 
the other hand r max (Ei) is very large, the great majority 
of accessible states will have higher energies, because the 
density of states usually decreases rapidly with energy. In 
consequence, there must be a value of r max (Ei) in-between 
that shows the desired property. In order to find this value 
we modify the radii after any proposed update v — t fi ac- 
cording to 



(l-e)r max (£y,if EvKEp 
(I + 2e)r max (£ ? ;),ii E v > 



(22) 



with Ei < E v < Ei+i and < e <C 1. It is easy to see 
that r max (i?i) will remain approximately unaltered if it 
has the desired characteristics, i.e., if E v < E^ in 66.6% 
of all cases. If the fraction of proposed moves leading to 
higher energies is too big, r max will be reduced and if it 
is too small, r max will be increased. In our simulation we 
used e = I0~ 3 . . . I0~ 5 and found little difference in perfor- 
mance. As expected, higher values of e allow faster conver- 
gence but lead to more noise in r mayi (Ei). However, in all 
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Figure 1: Maximal step length r max (i5) after a preliminary tun- 
ing procedure for N = 309 (E = -1820.684). 



Figure 2: Density of states for N = 309 (Eo = -1820.684), 
covering more then 3000 orders of magnitude. 



considered cases r max (Ei) converged quickly and caused 
update acceptance rates above 60% for all energies. In 
Fig. [TJ the obtained radii for the homopolymer of length 
N = 309 are shown. The used ratio 1:2 was chosen for the 
sake of simplicity. Different values might be found to be 
appropriate as well. The only restriction is that the de- 
sired fraction of updates to higher energies must be larger 
than 1/2. 

If the applied algorithm is able to find the valley of 
the global energy minimum, in principle the optimization 
allows us to come arbitrarily close to the ground state. Re- 
maining problems are of "technical" nature and consider 
the resolution of the energy scale and limits of numerical 
data types. In Fig. [2l the density of states g(E) for the 
309mer as obtained from two simulations is shown. After 
we investigated the general behavior and covered approx- 
imately 2000 orders of magnitude in the density of states, 
we resamplcd the region E < —1815 with a much higher 
energy resolution gaining further 1000 orders of magnitude 
mg{E). 



In a similar approach, attempted some time ago 14j . 
the authors applied analytic functions r max (efc) depending 
on the energy of the single particle k that is to be moved 
within a canonical ensemble. In contrast to the results pre- 
sented here, decisive improvements could not be achieved. 
Most likely this is in the first place due to the fact that in 
the canonical ensemble the potential for speedups is much 
smaller than in the multicanonical ensemble. We also be- 
lieve for two reasons that the energy of a single particle 
as the argument of r max is in general less favorable then 
the energy of the entire system. First, when the system 
approaches the ground state, the particles might possess 
differing energies but r max has to be close to zero for all of 
them. Secondly, the same displacement will cause smaller 
relative changes for the global energy than for the single- 
particle energy and, therefore, smaller changes in r max . 
Thus, the correction factor will be closer to unity if the 



global energy is used and the general acceptance will be 
higher and/or larger steps are possible. 

Notice that the described tuning procedure leads to a 
violation of the detailed balance condition which seemed to 
be of little relevance, though, presumably since e is small. 
Of course, the tuning must be ceased for the production 
run, in order to exclude this source of systematic error. 

If the considered system has continuous degrees of free- 
dom, this optimization procedure should in principle al- 
ways be applicable to basic Monte Carlo moves. However, 
the method might not work as described in the excep- 
tional situations when the density of states decreases with 
increasing energy. In these (rare) cases one should not rely 
on the proposed energy, but on the density of states itself, 
i.e., the radius has to be reduced (increased) if the update 
leads to an energy with a higher (lower) density of states. 
This was not necessary for the here investigated polymer 
model and since the density of states is not known a priori 
we employed the energy as reference. 

4-2. Bond-exchange moves 

While performing bond-exchange moves the positions 
of the monomers remain unchanged, but the bonds be- 
tween them are rearranged. In the past this type of confor- 
mational update has been applied mainly to lattice poly- 
mers [lj|, and applications for off- lattice polymers have 
also been documented and proven to be efficient 
For the sake of completeness we present the two different 
types used in our investigations. 

The first version, depicted in Fig. |3l consists of a swap 
of bonds between four nearby monomers. Initially, the 
monomers are labeled by numbers according to their posi- 
tion along the chain. Assuming two bonds have been cho- 
sen to be swapped, only one way exists to reconnect the 
chain without splitting the polymer. Let the contributing 
monomers be on the positions + and j + 1 with 
j > % + 1. It is obvious that if the iih bond between 
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Figure 3: Bond- exchange update. 



monomer i and i + 1 and the jth bond between monomer 
j and j + 1 are removed, different bonds can only be estab- 
lished between the ith and the jth monomer on the one 
side, and between the (i + l)th and the (j + l)th monomer 
on the other. Creating a bond between the (i + l)th and 
the jth monomer would result in a closed loop, since both 
monomers are already connected by a sequence of bonds. 
In our simulations, we first randomly choose an arbitrary 
bond i and determine afterwards which other bonds can 
possibly participate in an exchange update. Since in the 
employed model the bond length has an upper and a lower 
limit, only a few bonds are candidates. From this group 
the second bond j is then drawn randomly and the accep- 
tance probability is calculated. 

At this point it is important to recognize that also for 
this type of Monte Carlo move the probability for selecting 
the update, which is inversely proportional to the product 
of the number of bonds -/Vbonds an d the number of possible 
exchange partner bonds n" 1 *, often differs from that of the 
backward update. Both must be calculated and used for 
the determination of the acceptance probability according 
to (fTS|) . One obtains 



Wf, 



■ (, W 
mm 1, 



min 1, 



Pu ■ (AWrff)- 1 



= min 1, 



Pan 



cpb s 



P„n 



cpb 



(23) 



The order of monomers and bonds gets changed during 
the update and eventually appears to be totally random, 
if it is not restored by relabeling. 

If only the update just described is used, an end mono- 
mer will always remain an end monomer and the simula- 
tion would still be inefficient. Hence, we applied a second 




Figure 4: End-bond- exchange update. 



bond-exchange move shown in Fig. 2) Thereby we connect 
an end monomer to another nearby monomer and break 
the created loop by removing the old bond next to the 
formed junction. More explicitly, if we connect the first 
monomer to the jth, we obtain a ring of bonds connecting 
the first j monomers with a side chain branching off at the 
jth monomer. To remove the junction we have to delete 
the bond between the (j — l)th and the jth monomer. In 
the second case where the Nth monomer gets connected to 
the jth, the bond between the monomers j and j + 1 has 
to be deleted. Within the simulation we choose one of the 
end monomers and determine all monomers that are pos- 
sible partners for the update. Again, we draw monomer 
j from this set and, in order to be able to calculate the 
acceptance probability W a , it is necessary to consider the 
selection probabilities for the update in both directions: 



= min 1 



(24) 



with and n l 



being the numbers of possible ex- 
change partner monomers and i G {1, N}. 

The application of the two bond-exchange updates sig- 
nificantly increased the performance of the simulation and 
allowed larger changes of the polymer's configuration also 
in the "frozen" low-temperature regime. Even if there are 
no noticeable changes in monomer positions, the bonds are 
still quite flexible and arrange in a specific order when zero 
temperature is approached. Exemplified for the lowest- 
energy conformation of the 309mer, the length of each 
bond is shown in Fig. [SJ where the shell to which it belongs 
is represented by the symbol and the color. As a result of 
the icosahcdral packing, neighboring monomers are closest 
when they belong to neighboring shells. This makes these 
monomer pairs unfavorable for bonds, and in consequence 
only one bond each connects the inner shells, and at low 
T one end of the polymer is always located in the center. 

4-3. Monomer cut- and- paste update 

Below the liquid-solid transition the representative con- 
formations differ not only in the arrangement of the bonds, 
but in monomer positions as well. Even if the ground state 
is a perfect icosahedron, single monomers can be displaced 
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Figure 5: Bond ordering for the putative ground state of the 
309mer icosahedron. 



at low temperatures, thereby creating multiple surface de- 
fects (Fig. [6]). Transitions between these microstates can- 
not be performed with simple monomer displacements and 
bond-exchange moves only, since high energy barriers sep- 
arate favorable monomer positions on the surface of the 
icosahedron. Hence, we developed a fourth type of Monte 
Carlo move (Fig. [7]) to overcome this difficulty For this 
update, a monomer i is selected whose neighbors arc at a 
appropriate distance to be bonded themselves. In order to 
possess two neighbors the chosen monomer must not be an 
end of the polymer (1 ^ i ^ N). The position (r, (j>, z) s of 
monomer i is then determined according to a cylindrical 
coordinate system s defined as follows: The z-axis points 
through the neighboring monomers i — 1 and i + 1 and the 
origin is located in their midpoint. The further orienta- 
tion of s is irrelevant, because the original angle <f> will not 
be needed in the following. Now, monomer i is cut and a 
bond connecting the monomers i — 1 and i + 1 is created 
while another existing bond is removed in order to paste 
monomer i at its position. For that purpose, a second co- 
ordinate system s' is defined similar to s but based on the 
adjacent monomers of the removed bond, say monomer j 
and j + 1 (j 7^ i ^ j + 1). The coordinates r and z are 
now transposed from s to s' and a new angle <f> is drawn 
randomly from [0,27r). Again the angular orientation of 
s' can be arbitrary. Monomer i is now placed at this new 
position (r,<j)',z) s i and connected to the monomers j and 
3 + 1. 

The selection probabilities of the move and its inversion 
are identical, and no correction needs to be applied at this 
point. However, it is appropriate to introduce restrictions 
to the choice of the monomer to be moved and the bond 
to be split. If the polymer occupies a compact shape, the 
update has only a good chance of acceptance when per- 
formed at the surface, since moving a monomer within the 
interior, as well as from the center to the surface, implies 
a large increase in energy and a very low acceptance rate. 




Figure 7: Monomer cut-and-paste update. 



It is therefore useful to choose only bonds and monomers 
that are in regions of minor density, e.g., at the surface of 
a compact conformation. To estimate the density we use 
the number of contacts of a monomer (for details see 0j), 
i.e., the number of monomers to which its distance does 
not exceed a certain threshold. Since inner monomers at 
low temperature always have 12 contacts, we choose only 
monomers with less than 11 contacts and bonds that con- 
nect monomers with less than 12 neighbors. Unfortunately 
this leads to unequal selection probabilities and requires 
once more the introduction of a correction term. If n v m is 
the number of monomers to choose from and n\ i is the 
number of available bonds, we obtain 

w ^=^W¥\- (25) 

Here, n\ )i depends on i in a non-trivial way since bonds 
adjacent to monomer i must not be chosen. An alterna- 
tive way would be to allow choosing these bonds, but to 
immediately reject the update, once they are selected. 

5. Extensions to the Multicanonical Sampling Al- 
gorithm 

In the previous section we described how to overcome 
the problem of energy barriers through avoiding them by 
the application of certain update procedures, which is pos- 
sible in the described cases since the configurations on 
both "sides" of the barriers are rather similar. For the 
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bond-exchange update the monomer positions are identi- 
cal, and in the case of the cut-and-paste update, only a sin- 
gle monomer is moved. However, other barriers of different 
nature exist, and need to be treated with other strategies. 
As we have shown Q, the polymers adopt different ge- 
ometries corresponding to different optimizing strategics, 
resembling the behavior of atomic LJ clusters. This sim- 
ilarity has been already reported for a slightly different 
model [lH some time ago and is the result of the match- 
ing minimum distances of the two interaction potentials, 
which ensure that configurations minimizing the Lcnnard- 
Joncs potential also lead to low bond energies. Clusters 
and polymers both favor icosahedral crystal-like confor- 
mations at temperatures below the liquid-solid-transition. 
These conformations divide into two subgroups according 
to the type of the outer layer which can be either Mackay 
(fee) or anti-Mackay (hep) [l9| . Transitions between these 
two types occur at different temperatures, and for cer- 
tain system sizes, the investigation with standard Monte 
Carlo methods is difficult or impossible due to high free- 
energy barriers between different solid phases associated 
with Mackay or anti-Mackay growth. A second type of 
solid-solid transition that occurs for special system sizes in- 
volves non-icosahedral ground-state conformations, which 
can be of fee-, decahedral, or tetrahedral structure. These 
systems change to an icosahedral solid state at very low 
temperatures, posing a considerable challenge to the ap- 
plied simulation method. 

5.1. "Grand-multicanonicaV simulation 

First, we will consider the Mackay-anti-Mackay tran- 
sition within the surface of an icosahedral conformation. 
As already mentioned, the investigated LJ homopolymer 
behaves very similar to atomic LJ clusters. In the inter- 
val A £ [13,147], we find anti-Mackay ground states for 
13 < A < 31 and 55 < A < 81 while for the remaining 
polymer lengths Mackay ground states are favored. Ex- 
ceptions are A = 38, 75 - 77, 86, 87 d. 

Most of the systems with Mackay ground states un- 
dergo a transition to anti-Mackay conformations at a tran- 
sition temperature which generally increases with system 
size (Fig. [SJ. It turned out that this transition compli- 
cates the investigation, if it takes place at low tempera- 
tures, as for A = 31, or if the system is large, e.g., for 
A > 81. If standard methods like parallel tempering 0], 
multicanonical sampling Q or the Wang-Landau method 
[|| are applied, the system has to cross the barrier be- 
tween the Mackay and the anti-Mackay state many times 
in order to produce precise results. It turned out that 
this can be avoided by allowing the system to move also 
in A-direction, i.e., to change its size, during the simula- 
tion. The system is then able to circumvent the Mackay- 
anti-Mackay transition by changing A, and performing 
two liquid-solid transitions (Fig. [5$, which happens more 
frequently than the crossing of the Mackay-anti-Mackay 
transition line for sizes 81 < A < 110. 
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Figure 8: Sketch of the conformational state space at low temper- 
atures with "paths" to avoid the Mackay-anti-Mackay barrier. 



To move in A-direction we need a new Monte Carlo 
update that changes the system size at runtime. For- 
tunately, the monomer cut-and-paste update introduced 
above can be used as a starting point. If an increase of 
system size should be proposed, a bond k can be picked 
and coordinates of the new monomer are randomized. We 
again apply a cylindrical coordinate system s, defined by 
the adjacent monomers of the chosen bond: the z-axis 
points through these monomers and their midpoint de- 
fines the origin. The angular orientation is arbitrary. The 
coordinates (r,cj>,z) s have to be determined in order to 
be uniformly distributed in the hollow cylinder defined by 
?"min, r m ax, and Zmax (Fig-EJ). Therefore, <f> and z are drawn 
from constant distributions over the intervals [0, 2-k) and 



), respectively. Within 



' mm ; ' max 



) the de- 



sired probability density Pr{t) has to be proportional to 
the area of the cylinder shell with radius r, i.e., propor- 
tional to r itself. If we regard the radius r as a monotonic 
function of a uniformly distributed random number £: 

r = r(0, (26) 
where the probability density of £ is given by 

m) = {l : l^- (27, 

the fraction of points in the ring between r m ; n and r 
ir(r 2 — r 2 ■ ) 

v mm/ 



f(r) = 



ir(r 2 — r 2 ■ ) 

\ max mm/ 



has to equal £, 

e=/(r), 

which leads to 



V (''max '"min)£ 



(28) 



(29) 



(30) 



The inverse update meaning the reduction of the sys- 
tem size is simpler to accomplish. A monomer, which must 
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Figure 9: Coordinates for adding a new monomer. 




300 



250 



not be an end monomer, is chosen randomly and once more 
the coordinates (r' , <fi' , z') s i in a cylindrical system s' de- 
fined by the neighbors are determined. The update may 
only be performed if \z'\ < z max and r min < r' < r max , 
since otherwise the inverse update would be impossible, 
violating detailed balance. Note that in its present form 
the update contains another imbalance, since for the first 
choice the number of alternatives differs. If the system size 
should be increased, we choose from N — 1 bonds while, if 
the size is to be decreased, there are only N' — 2 monomers 
(with N' = N + I) to choose from. However, this imbal- 
ance can be neglected, since it does not effect the balance 
of conformations with identical N. 

To calculate the acceptance probability we first need 
the probability of each conformation. Again, we use a 
weight function uj{E, N) to produce a flat distribution but 
now in the two directions N and E. It is 



P w oc^£({X}),AT({X})) 
and with (|10p we easily obtain 



W. 



{X}^{X'} 



;(g({X'}),iV(X'))iy { 5 x , w{x} 
*(£({X}),iV(X))W? x 



{X}^{X'} 



(31) 



(32) 



Again, it is appropriate to choose only bonds and mono- 
mers from the surface. The adaptation of the method and 
the determination of W s are very similar to the procedure 
we discussed for the cut-and-paste update and are not re- 
peated here. Note that the imbalance mentioned in the 
last paragraph is cured this way, too. 

This algorithm proved to be surprisingly efficient. While 
it appeared to be impossible to investigate the full behavior 
of the lOOmer with standard multicanonical simulations, 
the simultaneous sampling of all chains with N < 147 
did not pose any major difficulties. Furthermore, we were 
able to derive the thermodynamics for all polymers of size 
13 < N < 309 down to T « 0.05 within a single simulation 
on a single Intel Xeon core (3.06GHz). This simulation in- 
volved 2 x 10 12 single updates and ran for approximately 
5 months. Some results are shown in Fig. 1101 

Note that this method is primarily not designed to in- 
vestigate the grand-canonical ensemble. Here, the focus is 
still on systems of fixed size and the merit lies in greater 
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Figure 10: Results from a single grand-multicanonical simula- 
tion: (a) specific heat, (b) temperature derivative of the nor- 
malized radius of gyration. 



efficiency in sampling them simultaneously and not in a 
physical understanding of polymerization processes. 

5.2. Multicanonical simulation with multiple weight func- 
tions 

The existence of non-icosahcdral ground states for a- 
tomic LJ clusters of certain sizes has been known for a 
long time, but the identification of these ground states is 
still regarded to be a major challenge to the applied algo- 
rithm. On the other hand, the investigation of the associ- 
ated solid-solid transitions is even more complicated, since 
the goal is not only to reach the ground-state conformation 
but also to maintain detailed balance and to measure the 
density of states very precisely. To the best of our knowl- 
edge there has been only one successful attempt to solve 
the problem for the 98-atom cluster [2(|, which involved 
the construction of an artificial energy landscape based on 
the prior knowledge of low-energy conformations. Here, 
we present an extension to the multicanonical approach 
which allows for investigating the solid-solid transitions of 
LJ polymers and clusters, but at the same time is general 
enough to be of use in other cases, too. 

In 0, Hj], we used the number of icosahedral cells to 
introduce a parameter v that indicates the geometrical 
state of the system: With high reliability we found v = 
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Figure 11: Decomposition of the microcanonical ensembles ac- 
cording to the different values of the order parameter v. 



for unstructured and for non-icosahcdral states, v = 1 for 
icosahedral states with Mackay overlayer, and v = 2 for 
icosahedral states with anti-Mackay overlayer. 

For the 98mer with large cutoff (r c = ho) the decom- 
positions of the "microcanonical" ensembles according to 
this parameter are depicted in Fig. QT] Since the differ- 
ent values of v belong to very different structures, lines 
between the different domains in Fig. QT] can only be pen- 
etrated in the high-energy regime. Hence, any algorithm 
producing these microcanonical distributions (e.g., sim- 
ulated tempering, parallel tempering, the multicanonical 
method or the Wang-Landau technique) is prevented from 
finding the tetrahedral ground-state conformation, since 
the probability to pass through the bottle neck belonging 
to v = at E « —500 is by far too small. The solution 
is to balance the probabilities of the three subensemblcs 
by introducing single weight functions for each value of v. 
Based on the multicanonical approach (|16[) , we use 



P { x} <x« v({X })(£({X})) 
to derive the acceptance probability 



(33) 



{X}^{X'} 



l, ^txog(CT))W7^ W{ x} ] 



j vW (E({X}))Wl 



X}^{X'} 



The remaining task is to tune the multiple weight function 
u> to allow each geometry to participate equally at any en- 
ergy and to enable the system to reach the energies where 
the solid-solid transition takes place. 

Results of applications of this algorithm are reported 
in detail in Ref. [8(- 

6. Conclusions 

In this paper, we described methods used to investigate 
the behavior of flexible homopolymers in much more detail 
and at much lower temperatures than it was previously 
possible. 



With the energy-dependent step length we introduced 
a novel general optimization scheme for basic Monte Carlo 
moves for systems with continuous degrees of freedom 
which allows constantly high acceptance rates everywhere 
in energy space. Applying this procedure in combination 
with multicanonical sampling we were able to estimate the 
density of states over several thousands of orders of mag- 
nitudes. 

We then described two bond-exchange moves and de- 
monstrated that these updates allow the reordering of poly- 
mer bonds without alteration of monomer positions. Sub- 
sequently, with the monomer-jump update we introduced 
a novel Monte Carlo move which increased the efficiency 
of the simulation further in two ways. First, the update 
allows the tunneling of energy barriers in the solid phase 
and second, it performs larger changes in the unstructured 
globular and the random coil phase. 

By enabling variations in system size at runtime we 
extended the multicanonical ensemble. This led to an ad- 
ditional gain in efficiency since the thus modified algorithm 
was able to circumvent certain energy barriers or to pene- 
trate them where they are low, i.e., at their "weak" points. 
As a result we obtained information over the entire state 
space over a large size interval from a single simulation. 

Finally, confronted with the problem of broken ergodic- 
ity and low-temperature solid-solid transitions, we devel- 
oped a second extension to the standard multicanonical 
technique. Due to the application of additional weight 
functions it is possible to retain ergodicity and to reach 
"hidden" ground states by circumventing the "blocking" 
states at intermediate temperatures. Although we yet have 
demonstrated the potential of this methods for hompoly- 
mcrs only, it is a general approach and, in combination 
with suitable order parameters, it might lead to substan- 
tial progress in the investigation of many other systems as 
well. 
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